Introduction
Compared with the well-known model plants Arabidopsis thaliana and rice, Brachypodium distachyon
(B. distachyon)
is more closely related to many important cereal crops such as wheat, barley,
oats, corn, and so on (Brkljacic et al. 2011; Catalan et al. 2015). Therefore B. distachyon is a currently recognized model organism for
Gramineae, and its biological characteristics, such as easy cultivation,
smaller habit, and short growth cycle (Scholthof et al. 2018), are conducive to
in-depth research on the development, evolutionary biology and abiotic stress
response of temperate grain crops.
Unlike linear RNA,
circular RNA (circRNA) is a special member of the
non-coding RNA family. Since circRNA forms a covalent
closed loop structure and does not have a polyadenylic
acid tail, it is resistant to RNase, which can effectively degrade linear RNA (Zhao et al. 2019).
The discovery of circRNA has been for decades, and it
was once considered a by-product of splicing errors (Sanger et al. 1976; Kos et al. 1986). However, with the
continuous development of high-throughput deep sequencing technology and
bioinformatics tools, a variety of circRNAs have been
discovered and identified in many organisms, such as humans (Jeck et al. 2013),
fruit flies (Westholm et al. 2014), mice (Fan et al. 2015), zebrafish (Shen et al. 2017),
Arabidopsis (Ye et al. 2015), rice (Lu et al. 2015)
and soybeans (Wang et al. 2020). These researches confirmed that circRNA is ubiquitous in eukaryotes. The recent research on
circRNA is mostly concentrated in mammals (Memczak et al. 2013),
whereas research in plants is still in its infancy. The whole genome identification
of circRNA in plants was first carried out in Arabidopsis thaliana and Oryza sativa. Ye et al. (2015, 2016)
identified 6,012 and 2,806 circRNAs from Arabidopsis
leaves and rice seedling root tissues, respectively. In addition to these two
well-known model plants, circRNA has also been
identified in other monocot and dicot species, such as soybean (Wang et al. 2020),
barley (Behrooz et al. 2016), and wheat (Xu et al. 2019). Similar to the
expression pattern of circRNA in animals (Westholm et al.
2014; Fan et al. 2015), circRNA in plants also exhibits tissue specificity and
responds to environmental stress (Lu et al. 2015; Behrooz et al. 2016). For example, rice
has 27 different expressions of circRNA under
conditions of sufficient phosphate and starvation (Ye et al. 2015). There
are 163 circRNAs in tomato showing a cold-responsive
expression pattern (Zuo et al. 2016). Increasing evidences indicate that circRNAs may play significant roles in a variety of
biological processes, such as miRNA binding, protein binding and
transcriptional regulation (Hansen et al. 2013; Memczak et al. 2013; Chen et al. 2016).
Drought is one of the
most serious adversity stresses affecting plant growth. With the global
warming, the shortage of water resources has become more and more serious,
which has directly led to the expansion of arid areas and the increase of aridity.
In the long-term evolutionary process, plants have formed their own defense
response mechanisms to deal with various stresses including water privation,
such as the accumulation of osmotic adjustment compounds (Sattar et al. 2020b),
stoma regulation systems (Sattar et al. 2020a), and signal
transduction pathways (Iqbal et al. 2017). However, the response of plants to drought
stress is the result of multi-angle and multi-level interaction, including not
only the transcription level and translation level, but also the regulation of
the post-transcriptional level (Chen et al. 2017b). As a result, in
addition to mRNA, there are also a large number of non-coding RNAs involved,
such as miRNA, lncRNA and circRNA. The researches
about the former two in drought stress have been widely reported (Ma et al. 2019;
Nadarajah and Kumar 2019; Feng et al. 2020),
whereas research on how circRNA functions in drought
stress response and regulation is rare. Therefore, this study adopted
high-throughput technology to sequence the transcriptome of the B. distachyon leaves under water shortage stress, and used
bioinformatics tools to analyze and identify the sequencing results, intending to
provide theoretical basis for further understanding and improvement of drought
response in Gramineous crops.
Materials and
Methods
Plant material,
growing environment and drought treatment
B.
distachyon Bd-21 was sprouted in a clean petri dish for 1 day and then
transferred to a 4°C incubator for 8 days for vernalization. After
vernalization, these seedlings were planted in a pot, placed in a light
incubator at 25°C with photoperiod 16 h light/ 8 h dark. When the seedlings
grew to the 7-leaf stage, they were subjected to drought treatment. The control
samples were watered once every other day, and the experimental group was
treated with water restriction in a drought treatment for 7 days. After the
treatment, the leaves of the seedlings were quickly frozen in liquid nitrogen
and stored in an ultra-low temperature refrigerator at -80°C for subsequent
experiment.
Library
construction, sequencing and circRNA identification
RNA
sequencing is a transcriptomics research method based on next-generation
sequencing technology. Each step of sample detection, library construction, and
sequencing is strictly controlled to ensure the output of high-quality data.
First, extract the total RNA of the biological sample, according to the Feng et al. (2020)
description. Then remove the rRNA and linear RNA. The RNA obtained is
theoretically only circular RNA, and then through the steps of reverse
transcription and PCR, and finally the enriched cDNA is subjected to
high-throughput sequencing. After base calling analysis, the original image
data files obtained by Illumina HiSeqTM2500 sequencing platforms are converted
into sequenced reads, namely raw reads. CircRNA from
the sequencing results could be identified through find circ (Memczak et al. 2013),
and the candidate circRNAs whose read count greater
than or equal to 2 were taken as the identified circRNA.
Expression
analysis of circRNAs
The
expression level of known and new circRNA in each
sample was counted, which was normalized by transcript per million tags (TPM)
(Zhou et
al. 2010). Differentially expressed genes sequencing (DEGseq) (Wang et al. 2010) was employed to
analyze the differences in circRNA expression between
different samples. The differential circRNAs were
screened from two aspects of fold change and corrected significance level (q
value). The default differential circRNA screening
condition is: q value < 0.01 and | log2 (fold change)| > 1.
Verification
of circRNA
Since
the splice junction mapping on the genome after reverse splice of circRNA is reversed, the specific primers for detecting circRNA are generally designed as divergent primers, which are
correspond to the convergent primers used to detect linear genes. Divergent
primers include two types. Type I is generally designed back-to-back so that
the splice junction is included in the middle of the product (Fig. 1A). Type II
is the junction overlapping divergent primers (JOD primers), which span the
splice junction, that is, the splice junction needs to be placed at the 3'end
of a primer so that the splice junction is included on the primer rather than
included in the PCR product (Fig. 1B). During PCR, JOD primers can only be
amplified normally when they match the target circRNA,
and it cannot be amplified when there is a mismatch, which greatly improves the
specificity. Convergent primers were also designed as positive controls for
linear transcripts (Fig. 1C). The designed primers were amplified and verified
using cDNA and gDNA as templates. In order to confirm the reliability of the
transcriptome sequencing results, 8 circRNAs were
randomly selected for qRT-PCR to verify their
expression patterns. The highly specific JOD primers were used as primers, and
the fluorescence quantification system and procedures were described by Feng et al. (2019)
and Ma et
al. (2019).
Bioinformatics
analysis of parent genes of circRNAs under drought
treatment
The
input data of circRNA differential expression was the
read count obtained in the analysis of circRNA
expression level. The software psRobot and psRNATarget were adopted to predict the miRNA binding site
of circRNA (Hansen et al. 2013) and the target gene
of miRNA (Dai and Zhao 2011),
respectively. After the differentially expressed circRNAs
between each group were acquired, according to the corresponding relationship
between circRNA and predicted miRNAs with circRNA binding sites, as well as target genes downstream
of miRNAs, the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes
(KEGG) enrichment analyses were carried out (Zhang et al. 2017). GO
(http://www.geneontology.org/) is an international standard classification
system for gene function. KEGG is the major public database about pathway (Kanehisa et al.
2008). Pathway significant enrichment analysis takes KEGG pathway as
the unit and applies hypergeometric test to find the pathway that is
significantly enriched in the target gene compared with the entire genome
background.
Results
Sequencing results statistics and circRNA
identification
To identify the circRNAs in B. distachyon
responsing to drought stress, we constructed the
RNA-Seq library from drought-treated B.
distachyon leaves and performed Illumina HiSeq high-throughput sequencing. The above library
construction, sequencing and bioinformatics analysis were carried out by
Beijing Novogene Bioinformatics Co., Ltd. Each step
of sample detection, library construction, and sequencing is strictly
controlled to ensure the output of high-quality data. The quality assessment of
the original data of the sample sequencing output is shown in Table 1. The CK
and clean reads under water loss are 122,429,492 and 111,457,314, respectively.
Clean bases reached 18.36G and 16.72G, respectively. The length distribution of
circRNA of all samples is shown in Fig. 2. It can be
seen from the Fig. 2 that the length of circRNA
obtained by sequencing is mainly concentrated below 1000nt, and the number
gradually decreases as the length increases. According to the source of circRNA in the genome and its constituent sequence, it can
be divided into exonic circRNAs,
intronic circRNAs and intergenic circRNAs.
Fig. 3 counts the sources of circRNA in CK and water
scarcity. Among them, There were 870 and 1,111 exonic
circRNAs in CK and drought, respectively, which are
dominated (53.28 and 52.85%), and followed by intronic circRNAs
with 710 and 920, and intergenic circRNAs are the
least with 53 and 71, respectively.
Analysis of
differentially expressed circRNAs under drought
The
expression level of known and new circRNA in each
sample was calculated, and the expression level was normalized by TPM (Zhou et al. 2010).
The differential circRNA was screened through the
fold change and the corrected significance level (padj/qvalue). Fig. 4 (volcano map) can infer the overall
distribution of different circRNAs. Among them, 332
and 265 were up-regulated and down-regulated separately, and 140 were not
dramatically different. The Venn diagram of differential expressed circRNAs under CK and drought treatment is shown in Fig. 5.
There are 229 and 303 circRNAs specifically expressed
under CK and drought, respectively, and 204 are co-expressed.
Validation of circRNA in B. distachyon
Fluorescence
quantitative polymerase chain reaction (qRT-PCR) was
performed to validate the 8 circRNAs randomly
selected from the circRNA-Seq analysis. Using highly
specific JOD primers to amplify circRNA, cDNA and
gDNA were used as amplification templates, and gDNA was used as negative
control (Fig. 6A). Since the JOD primer spans the circularization site, the
corresponding fragment cannot be amplified by using gDNA as a template. β-actin
was used as an internal reference gene to normalize the expression of these 8 circRNAs. The qRT-PCR results in
Fig. 6B showed that the expression levels of these 8 circRNAs
were consistent with the RNA-seq results, indicating the reliability of the RNA-seq
results. All the primers were showed in the Supplementary Table 2.
miRNA binding
site analysis
The
circRNA can inhibit the function of miRNA by binding
to miRNA (Hansen et al. 2013), that is, it can regulate gene expression by
acting as a miRNA sponge. Therefore, the miRNA binding site analysis of the
identified circRNA will help to further study the
function of circRNA. To investigate whether the circRNA in B. distachyon can affect the post-transcriptional
level of target genes by binding to miRNAs under dehydration treatment, a
bioinformatics method was employed to identify potential miRNAs target position
of circRNA that are differentially expressed in B. distachyon.
In
this study, 150 differentially expressed circRNAs have
putative miRNA binding sites, and each circRNA has at
least one predicted miRNA binding site, and a single miRNA can be targeted by
multiple circRNAs, and a single circRNA
can also target different miRNAs (Supplementary Table 1). For example, the
predicted bdi-miR5169a is targeted by 9 circRNAs,
while bdi_circ_0000061 can target 3 miRNAs in B. distachyon. In addition, based on predicted
miRNAs with binding sites, target genes downstream of miRNAs were screened. In
line with function annotations, they are divided into seven categories:
hormone-related, photosystem, stress response, regulation factor, ion
transport, transcription factor (TF), and ribosomal protein (Table 2,
Supplementary Table 2). Among the circRNAs
corresponding to these seven types of target genes, the number of up-regulated
expression was remarkably more than that of down-regulated expression (Table 3,
Supplementary Table 2).
%201355-1366_files/image002.png)
Fig. 1:
Schematic diagram of divergent and convergent primers. A, type I
divergent primer, the junction site is in the middle of the PCR product; B, type Ⅱ divergent
primer, overlapping with the junction site; C, convergent primer, used to verify linear transcripts. The black
rectangle represents the junction site
%201355-1366_files/image004.jpg)
Fig. 2: The
length distribution of circRNAs. circRNAs
below 1000 nt are 132 and 159 in CK and drought
samples, accounting for 40% and 43%, respectively
%201355-1366_files/image006.jpg)
Fig. 3: Statistical analysis of exonic circRNAs, intergenic circRNAs,
and intronic circRNAs in each sequenced sample. The
above three circRNAs are 870, 53, 710 in CK, and
1111, 71, 920 in drought sample, respectively
%201355-1366_files/image008.png)
Fig. 4: Volcano
plot of differentially expressed circRNA. The
abscissa represents the fold change of circRNA
expression in different samples, the ordinate represents the statistically
significant degree of circRNA expression change, the
scattered dots in the figure represent each circRNA,
and the blue dots represent circRNAs with no
significant differences, the red dots indicate significantly up-regulated circRNAs, and green dots indicate significantly
down-regulated circRNAs
%201355-1366_files/image010.jpg)
Fig. 5: Venn
diagram of differential expressed circRNA. The number
in each large circle represents the circRNA expressed
in this sample, and the overlapping part of the circle represents the circRNA co-expressed between samples
Bioinformatics
analysis of target genes predicted by circRNA under
drought
According to
the predicted miRNAs with circRNA binding sites,
target genes downstream of miRNAs were screened, and Gene Ontology and KEGG
enrichment analysis were performed on these target genes. The number of genes
in each GO term that were notably enriched is shown in a histogram (Fig. 7).
For biological processes, the category of metabolic process (GO: 0008152) is
the richest in GO terms. For cell components, the target genes of circRNA are mainly involved in cell (GO: 0005623) and cell
part (GO: 0044464). For molecular functions, the two most abundant categories are
binding (GO: 0005488) and catalytic activity (GO: 0003824). The scatter plot of
target gene KEGG
%201355-1366_files/image012.png)
Fig. 6:
Various experimental strategies validated the eight circRNAs
from transcriptome sequencing. A,
JOD primers (black back-to-back triangle pairs) and convergent primers (black
opposing triangle pairs) were adopted to amplify eight circRNAs
(bdi_0000061, bdi_0000429, bdi_0000041, bdi_0000461, bdi_0000672,
bdi_0000330, bdi_0000084, bdi_0000483) in cDNA and
gDNA. The former successfully amplified in cDNA but failed in gDNA, and the
latter worked on both cDNA and gDNA; B,
qRT-PCR validated the expression of the eight circRNAs
%201355-1366_files/image014.jpg)
Fig. 7: GO
enrichment histogram of target gene predicted by circRNA.
The abscissa is the GO term of the three major categories of GO, and the
ordinate is the number of target genes annotated to the term (including the
sub-terms of the term)
Table 1: Statistics of sequencing data
quality
|
Sample |
Raw reads |
Clean reads |
Clean bases |
Error rate (%) |
Q20 |
Q30 |
GC content (%) |
|
CK |
128,474,682 |
122,429,492 |
18.36G |
0.01 |
97.12 |
92.87 |
43.97 |
|
Drought |
116,766,272 |
111,457,314 |
16.72G |
0.01 |
97.21 |
92.92 |
44.70 |
%201355-1366_files/image016.png)
Fig. 8: KEGG
enrichment scatter plot of target gene predicted by circRNA.
The vertical axis represents the pathway, and the horizontal axis represents
the rich factor. The size of the dot indicates the number of target genes in
this pathway, and the color of the dot corresponds to different Q-value ranges
enrichment
is a graphical display of KEGG enrichment analysis results (Fig. 8). The degree
of KEGG enrichment is measured by rich factor, Q-value and the number of genes
enriched in this pathway. Among them, rich factor refers to the ratio of the
number of genes located in this pathway among differentially expressed genes to
the total number of genes located in this pathway among all annotated genes.
The greater the rich factor, the greater the degree of enrichment. Q-value is
the p-value after multiple hypothesis testing and correction. The range of
Q-value is [0, 1], which is closer to 0, the more significant. We selected the most
enriched 20 pathways for display in Fig. 8, among which biosynthesis of
secondary
metabolites, carbon metabolism and mRNA surveillance pathway are the three
pathways with the largest number of enriched genes.
Discussion
Previous studies believed that circRNA
was a byproduct of the transcription process (Cocquerelle et al. 1993; Memczak et al. 2013), but more and more
studies in the past decade have confirmed that circRNA
is a very important type of non-coding RNA, which plays a crucial role in a
variety of biological processes (Chen et al. 2017a). Most of the
researches on the function of circRNA have focused on
mammals (Salzman et al. 2012), however, circRNA
has also been identified in some plants, including Arabidopsis (Chen et al. 2017a),
rice (Lu et
al. 2015), wheat (Wang et al. 2017), and soybean (Wang et al. 2020).
In present study, transcriptome sequencing and identification of the circRNA in the leaves of B. distachyon under dehydrated stress were performed. A
total of 737 circRNAs were obtained, of which
Table 2: Function classification of the predicted
mRNA*
|
circRNA ID |
Binding site of mRNA |
target gene ID |
Function annotation |
|
|
Hormone related |
bdi_circ_0000045; bdi_circ_0000206; bdi_circ_0000214;
bdi_circ_0000257; bdi_circ_0000471; bdi_circ_0000472; bdi_circ_0000483;
bdi_circ_0000484; bdi_circ_0000485; bdi_circ_0000489; bdi_circ_0000559;
bdi_circ_0000576; bdi_circ_0000581; bdi_circ_0000582; bdi_circ_0000629;
bdi_circ_0000708 |
bdi-miR5174d-5p |
Bradi3g40830.1 |
auxin-responsive family protein |
|
bdi_circ_0000045; bdi_circ_0000047; bdi_circ_0000206;
bdi_circ_0000214; bdi_circ_0000257; bdi_circ_0000366; bdi_circ_0000471;
bdi_circ_0000472; bdi_circ_0000483; bdi_circ_0000484; bdi_circ_0000485;
bdi_circ_0000517; bdi_circ_0000559; bdi_circ_0000576; bdi_circ_0000581;
bdi_circ_0000582; bdi_circ_0000629; bdi_circ_0000708 |
bdi-miR5174f |
|||
|
bdi_circ_0000130; bdi_circ_0000131; bdi_circ_0000214;
bdi_circ_0000366; bdi_circ_0000471; bdi_circ_0000472; bdi_circ_0000576;
bdi_circ_0000581; bdi_circ_0000582 |
bdi-miR5174a |
|||
|
bdi-miR5174b-5p |
||||
|
bdi-miR5174e-5p.1 |
||||
|
bdi_circ_0000214; bdi_circ_0000471; bdi_circ_0000472;
bdi_circ_0000576; bdi_circ_0000581; bdi_circ_0000582 |
bdi-miR5174c-5p |
|||
|
bdi_circ_0000045 |
bdi-miR7757-3p.1 |
Bradi1g46060.1 |
abscisic acid responsive elements-binding factor 2 |
|
|
bdi_circ_0000144 |
bdi-miR7708a-3p |
Bradi1g46060.2 |
||
|
bdi_circ_0000061 |
bdi-miR171c-5p |
Bradi3g45880.9; Bradi3g45880.5; Bradi3g45880.10;
Bradi3g45880.3; Bradi3g45880.4; Bradi3g45880.6; Bradi3g45880.8 |
auxin response factor |
|
|
bdi_circ_0000491; bdi_circ_0000492 |
bdi-miR5166 |
Bradi3g04920.20 |
||
|
bdi_circ_0000324; bdi_circ_0000325 |
bdi-miR9486b |
Bradi2g44490.2; Bradi2g44490.1 |
auxin-like 1 protein |
|
|
bdi_circ_0000015 |
bdi-miR9491 |
|||
|
bdi_circ_0000082 |
bdi-miR164f |
Bradi1g22830.1 |
Auxin-responsive GH3 family protein |
|
|
bdi_circ_0000330; bdi_circ_0000331 |
bdi-miR528-3p |
Bradi2g46190.3 |
Transcriptional factor B3 family protein /
auxin-responsive factor AUX/IAA-related |
|
|
Photosystem |
bdi_circ_0000491;
bdi_circ_0000492 |
bdi-miR5166 |
Bradi1g24870.6; Bradi1g24870.2; Bradi1g24870.5;
Bradi1g24870.4; Bradi1g24870.3; Bradi1g24870.1 |
light harvesting complex photosystem II |
|
bdi_circ_0000366 |
bdi-miR9495 |
Bradi1g77340.1; Bradi1g77340.3 |
Mog1/PsbP/DUF1795-like
photosystem II reaction center PsbP family protein |
|
|
bdi_circ_0000041; bdi_circ_0000072; bdi_circ_0000395;
bdi_circ_0000429; bdi_circ_0000521; bdi_circ_0000531; bdi_circ_0000659;
bdi_circ_0000660; bdi_circ_0000708 |
bdi-miR5169a |
Bradi4g19720.1 |
photosystem II reaction center protein A |
|
|
bdi-miR5169b |
Bradi4g19720.1 |
|||
|
bdi_circ_0000483; bdi_circ_0000484; bdi_circ_0000485;
bdi_circ_0000632; bdi_circ_0000708 |
bdi-miR159a-5p |
Bradi1g25430.4 |
||
|
Stress response |
bdi_circ_0000294 |
bdi-miR5064 |
Bradi1g28090.2 |
chloroplastic drought-induced stress protein
of 32 kD |
|
CircRNA ID |
Binding site of mRNA |
Target gene ID |
Function annotation |
|
|
bdi_circ_0000009 |
bdi-miR5281a |
Bradi2g28040.2 |
Drought-responsive family protein |
|
|
bdi-miR5281b |
||||
|
bdi_circ_0000041; bdi_circ_0000092 |
bdi-miR7726a-5p |
Bradi3g34540.1; Bradi3g34540.5 |
early-responsive to dehydration stress protein (ERD4) |
|
|
bdi_circ_0000092 |
Bradi3g34540.4; Bradi3g34540.3 |
|||
|
bdi_circ_0000308 |
bdi-miR827-5p |
Bradi3g34540.6; Bradi3g34540.3; Bradi3g34540.4;
Bradi3g34540.5; Bradi3g34540.1 |
||
|
Ion transport |
bdi_circ_0000617 |
bdi-miR9496 |
Bradi1g45940.7; Bradi1g45940.8; Bradi1g45940.2 |
high-affinity K+ transporter 1 |
|
bdi_circ_0000072 |
bdi-miR437 |
Bradi5g26820.1; Bradi4g01430.1 |
K+ efflux antiporter 1 |
|
|
bdi_circ_0000009 |
bdi-miR5281a |
Bradi4g01430.1 |
K+ efflux antiporter 3 |
|
|
bdi-miR5281b |
||||
|
bdi_circ_0000617 |
bdi-miR9496 |
|||
|
bdi_circ_0000491; bdi_circ_0000492 |
bdi-miR5166 |
Bradi1g37860.9; Bradi1g37860.5; Bradi1g37860.6;
Bradi1g37860.7; Bradi1g37860.8; Bradi1g37860.11; Bradi1g37860.10;
Bradi1g37860.1; Bradi1g37860.12 |
K+ efflux antiporter 4 |
|
|
bdi_circ_0000708 |
bdi-miR169d |
Bradi1g76640.5 |
K+ efflux antiporter 5 |
|
|
Regulation factor |
bdi_circ_0000461 |
bdi-miR396e-5p |
Bradi1g09900.2; Bradi3g52547.1; Bradi5g20607.1;
Bradi1g50597.1; Bradi3g57267.1; Bradi3g39620.3; Bradi3g39630.1;
Bradi3g39590.1 |
growth-regulating factor |
|
|
bdi_circ_0000061 |
bdi-miR171c-5p |
Bradi2g23250.1 |
heat shock protein 70 |
|
|
bdi_circ_0000254 |
bdi-miR5184 |
Bradi5g24410.6; Bradi5g24410.3; Bradi5g24410.5;
Bradi5g24410.1; Bradi5g24410.2 |
jasmonate-zim-domain
protein 12 |
|
|
bdi_circ_0000444 |
bdi-miR5178-3p |
Bradi5g07080.3; Bradi5g07080.1; Bradi5g07080.4;
Bradi5g07080.2 |
Oxidoreductase, zinc-binding dehydrogenase family
protein |
Table 2: Continued
Table 2: Continued
|
|
CircRNA ID |
Binding site of mRNA |
target gene ID |
Function annotation |
|
bdi_circ_0000444 |
bdi-miR5178-3p |
Bradi5g07080.3; Bradi5g07080.1; Bradi5g07080.4;
Bradi5g07080.2 |
Oxidoreductase, zinc-binding dehydrogenase family
protein |
|
|
bdi_circ_0000461 |
bdi-miR7738-5p |
Bradi2g58130.5; Bradi2g58130.4; Bradi2g58130.6;
Bradi2g58130.3; Bradi2g58130.7 |
relative of early flowering 6 |
|
|
bdi_circ_0000228 |
bdi-miR7774-5p |
Bradi1g59560.4 |
splicing factor, putative |
|
|
bdi_circ_0000276 |
bdi-miR2118b |
Bradi5g22310.7; Bradi5g22310.10; Bradi5g22310.5;
Bradi5g22310.8; Bradi5g22310.4; Bradi5g22310.6 |
starch synthase 3 |
|
|
bdi_circ_0000367; bdi_circ_0000369; bdi_circ_0000370 |
bdi-miR7777-3p.1 |
Bradi5g22310.7; Bradi5g22310.10; Bradi5g22310.5;
Bradi5g22310.8; Bradi5g22310.4; Bradi5g22310.6; Bradi5g22310.3;
Bradi5g22310.9; Bradi5g22310.1; Bradi5g22310.2 |
||
|
TF |
bdi_circ_0000629 |
bdi-miR7782-3p |
Bradi3g08280.2; Bradi3g08280.1; Bradi3g08280.5;
Bradi3g08280.3; Bradi3g08280.4 |
basic helix-loop-helix (bHLH)
DNA-binding family protein |
|
bdi_circ_0000214; bdi_circ_0000471; bdi_circ_0000472;
bdi_circ_0000576; bdi_circ_0000581; bdi_circ_0000582 |
bdi-miR5174d-3p |
Bradi3g39927.1 |
||
|
bdi_circ_0000617 |
bdi-miR9496 |
|||
|
bdi_circ_0000366; bdi_circ_0000456; bdi_circ_0000489;
bdi_circ_0000490; bdi_circ_0000491; bdi_circ_0000511; bdi_circ_0000531;
bdi_circ_0000540 |
bdi-miR9493 |
Bradi1g54111.5 |
||
|
bdi_circ_0000072; bdi_circ_0000130; bdi_circ_0000131;
bdi_circ_0000242; bdi_circ_0000439; bdi_circ_0000511; bdi_circ_0000520;
bdi_circ_0000521; bdi_circ_0000607; bdi_circ_0000622; bdi_circ_0000623;
bdi_circ_0000624; bdi_circ_0000659; bdi_circ_0000660; bdi_circ_0000664;
bdi_circ_0000665; bdi_circ_0000682; bdi_circ_0000683; bdi_circ_0000733;
bdi_circ_0000734; bdi_circ_0000735 |
bdi-miR5171a |
Bradi1g05480.1 |
bZIP transcription factor family
protein |
|
|
bdi-miR5171b |
||||
|
bdi_circ_0000045 |
bdi-miR7757-3p.1 |
Bradi5g23340.3 |
||
|
bdi_circ_0000092 |
bdi-miR7743-3p |
Bradi1g46230.11; Bradi1g46230.8; Bradi1g46230.13;
Bradi1g46230.14; Bradi1g46230.12; Bradi1g46230.9 |
global transcription factor group A2 |
|
|
bdi_circ_0000041; bdi_circ_0000092 |
bdi-miR7726a-5p |
Bradi3g05260.1 |
K-box region and MADS-box transcription factor family
protein |
|
|
bdi_circ_0000708 |
bdi-miR169d |
Bradi3g54720.4 |
Plant-specific GATA-type zinc finger transcription
factor family protein |
|
|
bdi_circ_0000330; bdi_circ_0000331; bdi_circ_0000332;
bdi_circ_0000333; bdi_circ_0000334; bdi_circ_0000335 |
bdi-miR5198 |
|||
|
bdi_circ_0000396 |
bdi-miR5163a-3p |
Bradi1g08106.2 |
WRKY DNA-binding protein |
|
|
CircRNA ID |
Binding site of mRNA |
Target gene ID |
Function annotation |
|
|
bdi_circ_0000013; bdi_circ_0000041; bdi_circ_0000088;
bdi_circ_0000144; bdi_circ_0000162; bdi_circ_0000164; bdi_circ_0000165;
bdi_circ_0000166; bdi_circ_0000167; bdi_circ_0000211; bdi_circ_0000228;
bdi_circ_0000334; bdi_circ_0000335; bdi_circ_0000344; bdi_circ_0000366;
bdi_circ_0000458; bdi_circ_0000581; bdi_circ_0000582; bdi_circ_0000687;
bdi_circ_0000721; bdi_circ_0000736 |
bdi-miR5175b |
Bradi1g08106.4 |
||
|
bdi_circ_0000045 |
bdi-miR9494 |
Bradi2g19070.1 |
||
|
bdi_circ_0000061 |
bdi-miR171c-5p |
Bradi4g33370.2; Bradi4g33370.9 |
||
|
bdi_circ_0000330; bdi_circ_0000331; bdi_circ_0000672 |
bdi-miR528-5p |
Bradi4g28280.1 |
||
|
bdi_circ_0000045 |
bdi-miR9494 |
Bradi2g42023.1 |
||
|
bdi_circ_0000288; bdi_circ_0000366 |
bdi-miR164c-3p |
Bradi2g30695.1 |
||
|
bdi_circ_0000172 |
bdi-miR7771-3p |
Bradi2g30800.2 |
||
|
Ribosomal protein |
bdi_circ_0000172 |
bdi-miR7771-3p |
Bradi3g09780.1 |
Late embryogenesis abundant (LEA) hydroxyproline-rich
glycoprotein family |
|
bdi_circ_0000172 |
bdi-miR7771-3p |
Bradi4g34750.1 |
Ribosomal protein L7Ae/L30e/S12e/Gadd45 family protein |
|
|
bdi_circ_0000045 |
bdi-miR7757-3p.1 |
Bradi1g71200.1 |
Ribosomal protein S3Ae |
*Part of the classification results is
displayed. See Supplementary Table 2 for the complete results
exon circRNA exceeded 50%
(Fig. 3). The high ratio of exon circRNA is similar
with the Arabidopsis exon circRNA under high
temperature (Pan et al. 2018). It is speculated that the increase in the
number of exon circRNA may be due to adversity
stress. The increase in exon circRNA not only
responds to plant stress, but may also enhance stress resistance. It has been
confirmed that circRNA in plants can respond to a
variety of stresses. Zuo et al. (2016) verified that 163 circRNAs
in tomato respond to cold stress, while 62 circRNAs
in wheat are differentially expressed under drought stress compared to the
control (Wang
et al. 2017). There are 597 circRNAs with
significant differential expression in the leaves of B. distachyon
under Table 3: Statistics of circRNA
expression in seven functional categories
|
Function
classification |
circRNA |
Up-regulated circRNA |
Down-regulated
circRNA |
|
Hormone
related |
124 |
78 |
46 |
|
Photosystem |
47 |
27 |
20 |
|
Stress
response |
30 |
20 |
10 |
|
Ion transport |
34 |
21 |
13 |
|
Regulation
factor |
94 |
56 |
38 |
|
TF |
141 |
88 |
53 |
|
Ribosomal
protein |
130 |
79 |
51 |
dehydration stress, which is consistent with previous
research results, that is, circRNA can respond to
multiple stresses and may further participate in the regulation of stress
resistance.
The predicted target
genes of circRNAs were analyzed by bioinformatics.
The GO results showed that they were related to multiple functions, including
different biological processes, cellular components, and molecular functions. EGG
analysis enriched 20 pathways, most of which were involved in drought response.
For example, the pathway with the largest number of enriched genes is
biosynthesis of secondary metabolites, and there are as many as 89
hormone-related genes in the functional classification of target genes,
including auxin-responsive family protein, auxin-responsive factor, and
abscisic acid responsive elements-binding factor. Phytohormones are active
substances induced by plant cells receiving specific environmental signals.
They not only regulate plant growth and development, but also participate in a
variety of stress responses, including low temperature (Xin et al. 2019), water
deficit (Hu
et al. 2018), salt (Yu et al. 2020) and so on (Wang et al. 2018).
The exogenous application of ABA can increase the cold tolerance of rice plants
during the flowering period and thus increase the seed setting rate (Xiang et al. 2019).
Besides, exogenous application of gibberellin could relieve the adverse effects
of salinity, drought, and heat stresses on the growth of Capsicum annuum L., including increase biomass and chlorophyll
content, as well as improve photosynthetic efficiency (Khan et al. 2015). Since
circRNA can act as a molecular sponge of miRNA and
prevent it from regulating target mRNA and controlling gene expression (Memczak et al. 2013),
the expression trend of most circRNAs is positively
correlated with their corresponding target genes. Under drought treatment,
there were more markedly up-regulated circRNA in the
leaves of B. distachyon than down-regulated
expression. For example, the circRNAs predicted to be
hormone-related genes, 78 were up-regulated and 46 were down-regulated.
Therefore, it is speculated that the circRNA
corresponding to the hormone-related gene is involved in dehydrated response of
B. distachyon in a positive or negative
manner, and contributes to the regulation of its stress response.
In addition to
hormone-related genes, the circRNA corresponding to
target genes whose function is annotated as stress response, regulation factor,
photosystem and ion transport are also up-regulated more than down-regulated.
Drought is one of the common environmental stresses. Genes that respond to
drought stress have always been research hotspots. With the continuous development
of high-throughput sequencing technology, the function of non-coding RNA is
gradually recognized by people. Apart from miRNA and lncRNA, circRNA in plants also play a significant role in stress
resistance. The circRNA and its target genes obtained
in this analysis can be divided into seven categories: hormone-related,
photosystem, stress response, regulation factor, ion transport, TF and
ribosomal protein according to their functional classification. Among them,
photosystem- and hormone-pathway, as known regulation networks, participate in
plant drought response (Chen et al. 2016; Rao and Chaitanya 2016), including a series of
responses, like light harvesting complex photosystem II, photosystem II reaction
center protein A, auxin response factor and auxin transport protein genes, are
involved in dehydration response, which is consistent with previous studies (Chen et al. 2016;
Rao and Chaitanya 2016). Furthermore, Rai et al. (2012) confirmed that
early-responsive to dehydration stress protein (ERD4) is highly expressed in
plant drought response and could effectively improve the desiccation
adaptability of plants. ERD4 in the stress response classification was
up-regulated under drought treatment, which was in line with the results of Rai
et al. Under adversity stress, a large number of genes are up-regulated and
expressed to positively regulate the stress resistance of plants, whereas these
emergency response genes are inhibited by related mechanisms in normal growth
conditions. During the evolution process, plants have established related
mechanisms to inhibit the expression of stress-related genes under normal
conditions. The main role of growth-regulating factors (GRFs) is to inhibit the
stress response under non-stress conditions. Compared with wild-type
Arabidopsis and AtGRF7 overexpression
lines, the Arabidopsis atgrf7 mutant
is more tolerant to salt and dehydrated stress (Kim et al. 2012), therefore, GRF is
down-regulated under stress conditions to activate related stress genes to improve
plant resistance. The circRNAs, whose target genes
belong to the growth-regulating factor classification, were distinctly
down-regulated under drought treatment, which is consistent with the negative
regulation of plant stress resistance by GRF, which may be a crucial regulatory
network of stress response in the plant.
The K+
efflux antiporter (KEA), for K+, is functioning in maintaining
pH, the active accumulation and balance of K+ and ion homeostasis (Zhu et al. 2018).
Additionally, soybean GmKEAs
gene family up-regulated in the beginning and down-regulated later under excessive
potassium stress (Tao et al. 2015). The response of GmKEAs gene to abiotic stress
indicates that it is involved in soybean resistance. The target genes, whose
function annotated as KEA1, KEA3, KEA4 and KEA5 in B. distachyon,
their corresponding circRNAs were also up-regulated
after water shortage, which may also be involved in the response to osmotic
stress. WRKY TF exerts a critical role in a variety of biological processes, including
plant growth, development and adversity stress. Gao et al. (2018) overexpressed TaWRKY2 in wheat to enhance the dehydrated
resistance of transgenic plants. Shirazi et al. (2019) cloned the stress
response BnWRKY57 in Brassica napus and analyzed its
expression. Those results showed that BnWRKY57
not only responded to drought and salt stress, but also enhanced the resistance
of Brassica napus to the two
stresses. The circRNA whose target gene function was
annotated as WRKY DNA-binding protein was obviously up-regulated after drought
treatment, and it was speculated that it not only participated in the drought
response of B. distachyon, but also exerted an
influence in stress defense.
Moreover, some of
the other predicted target genes in the functional classification may also be
involved in the drought response of B. distachyon.
However, more evidence is needed to further prove its specific function, and
the mechanism of action may depend on circRNA
regulation.
Conclusion
Through transcriptome sequencing of B. distachyon leaves under drought
treatment, a total of 737 circRNAs were obtained in
this study, of which exon circRNA exceeded 50%.
Differential expression analysis showed that there were 332 and 265
up-regulated and down-regulated expressions, respectively and 140 had no
significant difference. There were 229 and 303 circRNAs
specifically expressed under CK and drought, respectively, and 204 were
co-expressed. Bioinformatics methods were used to identify potential miRNA targets
of differentially expressed circRNAs. Among them, 150
differentially expressed circRNAs had putative miRNA
binding sites. Based on the predicted miRNA with binding site, the target gene
downstream of the miRNA was screened. Analysis of GO and KEGG shows that the
categories of metabolic process (GO: 0008152), cell (GO: 0005623) and cell part
(GO: 0044464), binding (GO: 0005488) and catalytic activity (GO: 0003824) are
the most abundant entry in GO. Among the 20 most significant pathway entries,
biosynthesis of secondary metabolites, carbon metabolism and mRNA surveillance
pathway are the three pathways with the largest number of enriched genes.
According to the functional annotation of target genes, they are divided into
seven categories: hormone-related, photosystem, stress response, regulation
factor, ion transport, TF and ribosomal protein. CircRNA
is an important type of non-coding RNA in organisms, and these circRNAs involved in desiccation response may exert a
crucial role in drought response and defense of B. distachyon. The valuable information
also provides a theoretical basis for further identification of candidate genes
participated in drought regulation.
Acknowledgements
This research was financially supported by the Natural
Science Foundation of Henan Province (Grant No. 202300410133), the Key Research
Projects of Higher Education in Henan Province (Grant No. 19B210002), and Student
Research Training Program of Henan University of Science and Technology (Grant
No. 2020338 and 2020335).
Author Contributions
Yalan Feng performed the experiments and wrote the manuscript; Shuang Zhou and
Jun Zhang analyzed the experiment data and wrote the manuscript; Ke Xv, Yating
Li and Mengzhen Zhang sampled material and analyzed
the experiment data; Chao Ma and Youjun Li designed
the research and revised the manuscript.
Conflict of
Interest
Authors declare no conflict of interest.
Data Availability
The data and supplementary material are all available
online.
Ethics Approval
This research does not involve
the ethical approval.
References
Dai XB, XP Zhao (2011). psRNATarget: A plant small RNA
target analysis server. Nucl Acids Res
39:155‒159
Nadarajah K, IS Kumar (2019). Drought response
in rice: The miRNA story. Intl J Mol Sci
20:3766–3785
Scholthof KBG, S Irigoyen, P Catalan,
KK Mandadi (2018). Brachypodium: a
monocot grass model genus for plant biology. Plant Cell 30:1673‒1694
Shen YD, XW Guo, WM Wang (2017).
Identification and characterization of circular RNAs in zebrafish. FEBS Lett 591:213‒220
Wang XL, RR Qin, RH Sun, JJ Wang, XG Hou, L Qi, J Shi,
XL Li, YF Zhang, PH Dong (2018). No post-drought compensatory growth of corns
with root cutting based on cytokinin induced by roots. Agric Water Manage 205:9‒20
Ye CY, L Chen, C Liu, QH Zhu, L Fan (2015). Widespread
noncoding circular RNAs in plants. New
Phytol 208:88‒95